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10 ABSTRACT 

1 1 Mapping hybrid defects in contact zones between incipient species can identify genomic regions 

12 contributing to reproductive isolation and reveal genetic mechanisms of speciation. The house 

13 mouse features a rare combination of sophisticated genetic tools and natural hybrid zones 

14 between subspecies. Male hybrids often show reduced fertility, a common reproductive barrier 

15 between incipient species. Laboratory crosses have identified sterility loci, but each encompasses 

16 hundreds of genes. We map genetic determinants of testis weight and testis gene expression 

17 using offspring of mice captured in a hybrid zone between M. musculus musculus and M. m. 

18 domesticus. Many generations of admixture enables high-resolution mapping of loci contributing 

19 to these sterility-related phenotypes. We identify complex interactions among sterility loci, 

20 suggesting multiple, non-independent genetic incompatibilities contribute to barriers to gene 

21 flow in the hybrid zone. 
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22 INTRODUCTION 

23 New species arise when reproductive barriers form, preventing gene flow between 

24 populations (Coyne and Orr 2004). Recently, two approaches have substantially advanced 

25 understanding of the genetic mechanisms underlying reproductive isolation (reviewed in Noor 

26 and Feder 2006; reviewed in Wolf et al. 2010). Genetic crosses in the laboratory involving model 

27 organisms have identified loci and genes causing hybrid defects, a common type of reproductive 

28 barrier caused by genetic interactions between divergent alleles (Bateson 1909; Dobzhansky 

29 1937; Muller 1942). In nature, recent technological advances enable fine-scale characterization 

30 of genome -wide patterns of divergence between incipient species and variation in hybrid zones. 

3 1 For example, 'islands of divergence' have been reported in species pairs from 

32 taxonomically diverse groups (Turner et al. 2005; Nadeau et al. 201 1; Nosil et al. 2012; Ellegren 

33 et al. 2013; Hemmer-Hansen et al. 2013; Renaut et al. 2013; Carneiro et al. 2014; Poelstra et al. 

34 2014; Schumer et al. 2014). These high-divergence genomic outlier regions are sometimes 

35 referred to as 'islands of speciation,' resistant to introgression because they harbor genes causing 

36 reproductive isolation. However, other forces can create similar genomic patterns, thus islands 

37 may not always represent targets of selection that contributed to speciation (Noor and Bennett 

38 2009; Turner and Hahn 2010; Cruickshank and Hahn 2014). 

39 An alternative approach to identify genomic regions contributing to reproductive 

40 isolation is to map known reproductive barrier traits in naturally hybridizing populations. The 

41 potential for mapping in hybrid zones is long-recognized (Kocher and Sage 1986; Szymura and 

42 Barton 1991) (Harrison 1990; Briscoe et al. 1994; reviewed in Rieseberg and Buerkle 2002). 

43 Hybrid zones are "natural laboratories for evolutionary studies" (Hewitt 1988) enabling 

44 investigation of speciation in progress. The Dobzhansky-Muller model predicts that hybrid 

45 incompatibilities between incipient species accumulate faster than linearly with time (Orr 1995), 

p. 2 
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46 thus investigating taxa in the early stages of speciation facilitates identification of 

47 incompatibilities that initially caused reproductive isolation vs. incompatibilities that arose after 

48 isolation was complete. 

49 Despite these advantages, few studies have mapped barrier traits or other fitness-related 

50 traits in nature, due to the logistical challenges of collecting dense genome-wide genetic markers 

51 in species with admixed populations and well-characterized phenotypes. Examples include 

52 associations between pollen sterility and genomic regions showing reduced introgression in a 

53 sunflower hybrid zone (Rieseberg et al. 1999) and loci contributing to variation in male nuptial 

54 color and body shape mapped in a recently admixed stickleback population (Malek et al. 2012). 

55 House mice (Mus musculus) are a promising model system for genetic mapping in natural 

56 populations (Laurie et al. 2007), and have an abundance of genetic tools available to ultimately 

57 isolate and characterize the causative genes underlying candidate loci. Three house mouse 

58 subspecies - M. m. musculus, M. m. domesticus andM m. castaneus - diverged -500,000 years 

59 ago from a common ancestor (reviewed in Boursot et al. 1993; Salcedo et al. 2007; Geraldes et 

60 al. 2008). M. m. musculus and M. m. domesticus (hereafter musculus and domesticus) colonized 

61 Europe through different geographic routes and meet in a narrow secondary contact zone running 

62 through central Europe from Bulgaria to Denmark (Sage et al. 1986; Boursot et al. 1993). 

63 Genome -wide analyses of patterns of gene-flow in several geographically distinct transects 

64 across the hybrid zone have identified genomic regions showing reduced introgression, which 

65 may contribute to reproductive isolation (Tucker et al. 1992; Macho lan et al. 2007; Teeter et al. 

66 2008; 2010; Janousek et al. 2012). 

67 Reduced male fertility is common in wild-caught hybrids (Turner et al. 2012; 

68 Albrechtova et al. 2012) and in musculus - domesticus hybrids generated in the laboratory 

p. 3 
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69 (Britton-Davidian et al. 2005; reviewed in Good et al. 2008a), implying hybrid sterility is an 

70 important barrier to gene flow in house mice. Mapping studies using Fi, F 2 and backcross hybrids 

71 generated from laboratory crosses between house mouse subspecies have identified many loci 

72 and genetic interactions contributing to sterility phenotypes 

73 (Storchova et al. 2004; Good et al. 2008b; White et al. 201 1; Dzur-Gejdosova et al. 2012; Turner 

74 et al. 2014). Prdm9, a histone methyltransferase, was recently identified as a gene causing Fi 

75 hybrid sterility, and is the first hybrid incompatibility gene identified in mammals (Mihola et al. 

76 2009). Comparisons between different Fi crosses show that hybrid sterility alleles are 

77 polymorphic within subspecies (Britton-Davidian et al. 2005; Good et al. 2008a). Furthermore, 

78 reduced fertility phenotypes observed in nature vary in severity; complete sterility, as 

79 documented in some Fi crosses, appears to be rare or absent in the hybrid zone (Turner et al. 

80 2012; Albrechtova et al. 2012). The relevance of genetic studies of sterility in early-generation 

81 hybrids in the laboratory to understanding barriers to gene flow between later-generation hybrids 

82 in nature has yet to be verified (Britton-Davidian et al. 2005). 

83 Here, we map sterility-related phenotypes in hybrid zone mice to investigate the genetic 

84 architecture of reproductive isolation between incipient species. We performed a genome-wide 

85 association study (GWAS) to map testis weight and testis gene expression in 185 first generation 

86 lab-bred offspring of wild-caught hybrid mice (Figure 1-figure supplement 1). GWAS have been 

87 powerful in humans, loci contributing to hundreds of quantitative traits associated with disease 

88 and other phenotypic variation have been identified (reviewed in Stranger et al. 201 1). Examples 

89 of GWAS for fitness-related traits in non-humans are only beginning to emerge (Johnston et al. 

90 201 1; Filiault and Maloof 2012; Magwire et al. 2012). 
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91 Our hybrid zone GWAS identified genomic regions associated with variation in relative 

92 testis weight (testis weight/body weight) and genome-wide testis expression pattern, including 

93 regions previously implicated in hybrid sterility as well as novel loci. Motivated by the 

94 Dobzhansky-Muller genetic model of hybrid defects, we tested for genetic interactions 

95 (Dobzhansky-Muller interactions - "DMIs") between loci affecting testis weight or expression 

96 pattern. All loci except one showed evidence for interaction with at least one partner locus, and 

97 most interact with more than one partner. The deviations in phenotype associated with most 

98 interactions were large - affected individuals have phenotypes below the range observed in pure 

99 subspecies - suggesting these interactions indeed are hybrid incompatibilities. High resolution 

100 provided by mapping in natural hybrids enabled identification of one or few potential causative 

101 genes for many loci. 
102 

103 RESULTS 

104 Sterility-associated phenotypes. We investigated two phenotypes in males from the 

105 house mouse hybrid zone: relative testis weight (testis weight/ body weight) and genome-wide 

106 testis gene expression pattern. Both of these phenotypes have previously been linked to hybrid 

107 male sterility in studies of mice from crosses between musculus and domesticus and mice from 

108 the hybrid zone (Britton-Davidian et al. 2005; Rottscheidt and Harr 2007; reviewed in Good et 

109 al. 2008a; 2010; White et al. 201 1; Turner et al. 2012; 2014). We refer to these as "sterility 

110 phenotypes," following conventional terminology in the field, however, it is important to note 

1 1 1 that the severity of defects observed in most hybrid zone mice are consistent with reduced 

112 fertility/partial sterility (Turner et al. 2012). 

1 13 Testis expression PCI (explaining 14.6% of the variance) is significantly correlated with 

1 14 relative testis weight (cor = 0.67, P = 2xl0" 16 ) indicating there is a strong association between 

p. 5 
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1 15 those two sterility phenotypes (Figure 1 - figure supplement 2). Principal component 2 (PC2, 

116 8.1% variance) is strongly correlated with hybrid index (% musculus autosomal SNPs: cor = 

1 17 0.75, P = 2xl0" 16 ), thus the effect of hybrid defects do not obscure subspecies differences in 

118 expression. 

119 In the mapping population, 19/185 (10.2%) individuals had relative testis weight below 

120 the minimum observed in pure subspecies males and 21/179 (1 1.7%) individuals had expression 

121 PCI scores below (PCI = -46.97) the pure subspecies range. 

122 Association mapping. We identified 55 SNPs significantly associated with relative testis 

123 weight (false discovery rate (FDR)<0.1; Figure 1A), clustered in 12 genomic regions (of size 1 

124 bp - 13.3 Mb; Table 1). Three regions on the X chromosome were significant using a more 

125 stringent threshold determined by permutation. We report GWAS regions defined using a 

126 permissive significance threshold because we plan to combine mapping results from multiple 

127 phenotypes to identify candidate sterility loci, based on the idea that spurious associations are 

128 unlikely to be shared among phenotypes. Significant regions were located on the X chromosome 

129 and 9 autosomes, suggesting a minimum of 10 loci contribute to variation in testis weight. It is 

130 difficult to estimate the precise number of genes involved, because the extent of linkage 

131 disequilibrium (LD) around a causative mutation depends on the phenotypic effect size, 

132 recombination rate, allele frequency, and local population structure. Multiple significant regions 

133 might be linked to a single causative mutation, or conversely, a significant region might be 

134 linked to multiple causative mutations in the same gene or in multiple genes. 

135 We identified 50 genomic regions significantly associated with expression PCI, 

136 comprising 453 significant SNPs (Table 2, Figure IB) located on 18 autosomes and the X 
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137 chromosome. Five regions on the X and chromosome 1 were significant using a more stringent 

138 permutation-based threshold. 

139 To gain further insight into associations between sterility and gene expression, we 

140 mapped expression levels of individual transcripts. A total of 18,992/36,323 probes showed 

141 significant associations with at least on SNP. We focused on trans associations (SNP is located 

142 on different chromosome from transcript), based on evidence from a study in F 2 hybrids that 

143 trans expression QTL (eQTL) are associated with sterility while cis eQTL are predominantly 

144 associated with subspecies differences (Turner et al. 2014). To identify SNPs significantly 

145 enriched for trans associations with expression, we used a threshold set to the 95% percentile 

146 counts of significantly associated probes across all SNPs (30 probes, Figure 1C). 

147 There was substantial overlap between mapping results for testis weight and expression 

148 PC 1 ; 48/55 SNPs significant for relative testis weight (9 regions) were also significant for 

149 expression PCI . A permutation test, performed by randomly shuffling the positions of GWAS 

150 regions in the genome, provides strong evidence that this overlap is non-random (PO.0001, 

151 10,000 permutations). Most SNPs significant for testis weight and/or expression PCI were 

152 significantly enriched for trans associations with individual transcripts (relative testis weight: 

153 49/55 SNPs, 8/12 regions; PCI : 440/453 SNPs, 50/50 regions). The combined mapping results 

154 provide multiple lines of evidence for contributions of all 50 PCI regions and 9/12 testis weight 

155 regions. The three testis-weight regions (RTW04, RTW05, RTW08) not significantly associated 

156 with testis expression phenotypes are more likely to be spurious and are weaker candidates for 

157 future study. 

158 Genetic interactions. Power to identify pairwise epistasis in GWAS for quantitative 

159 traits is limited even with very large sample sizes, due to multiple testing issues (e.g. Marchini et 



p. 7 



Downloaded from http://biorxiv.org/on September 18, 2014 

Turner & Harr 

160 al. 2005). The Dobzhansky-Muller predicts that the effect of each hybrid defect gene depends on 

161 interaction with at least one partner locus. Hence, for hybrid sterility traits, there is a hypothesis- 

162 driven framework in which to limit tests for epistasis to a small subset of possible interactions. 

163 We tested for genetic interactions between all pairs of significant SNPs (FDR <0. 1) 

164 located on different chromosomes for testis weight and for expression PCI . We identified 142 

165 significant pairwise interactions for relative testis weight, representing 22 pairs of GWAS 

166 regions (Figure 2A). These results provide evidence for a minimum of 13 autosomal-autosomal 

167 and five X - autosomal interactions affecting testis weight. 

168 We identified 44,145 significant interactions between SNPs for expression PCI. The 913 

169 GWAS region pairs provide evidence that at least 144 autosomal-autosomal interactions and 18 

170 X-autosomal interactions contribute to expression PCI (Figure 2B). 

171 Effect size. We used deviations from population means for single SNPs and two-locus 

172 genotypes to estimate the phenotypic effects of GWAS regions and interactions (Figure 3A,B). 

173 As expected, interactions had greater effects, on average, than single loci for both pheno types 

174 (relative testis weight: single locus mean = -1.81 mg/g, interaction mean = -4.07 mg/g; 

175 expression PCI: single locus mean = -81.51, interaction mean = -130.77). We provide examples 

176 of autosomal- autosomal and X-autosomal SNP pairs with significant interactions for each 

177 pheno type in Figure 3 C. It is important to note that mean deviations are rough estimates of effect 

178 sizes, which don't account for family structure. 

179 It is possible that some of the GWAS regions we mapped contribute to quantitative 

180 variation within/between subspecies, rather than hybrid defects. The lowest genotypic means for 

181 most interactions fell below the range observed in pure subspecies (relative testis weight: 19/22 

p. 8 
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182 (86.3%) region pairs; expression PCI : 877/913 (96%) region pairs; Figure 3AB), consistent with 

183 the hypothesis that interactions represent Dobzhansky-Muller incompatibilities. 

184 Mapping simulations. We performed simulations to assess the performance of the 

185 mapping procedure for different genetic architectures by estimating the power to detect causative 

186 loci and the false positive rate. We simulated phenotypes based on two-locus genotypes from the 

187 SNP dataset using genetic models for nine genetic architecture classes (i.e. autosomal vs. X 

188 linked, varied dominance) with parameters based on the observed distribution of relative testis 

189 weight (Figure 4-figure supplements 1-2). 

190 The distribution of distances to the causal SNP for all significant SNPs located on the 

191 same chromosome (Figure 4-figure supplement 3) shows that the majority of significant SNPs 

192 (62.7%), are within 10 Mb of the causal SNP, however a small proportion of significant SNPs 

193 are >50 Mb from the causal SNP. In most cases, causal SNPs detected at long distances also had 

194 significant SNPs nearby, for example 83.4% of loci with significant SNPs 1-10 Mb distant also 

195 has significant SNPs within 1 Mb. These results provide support for our choice to define 

196 significant GWAS regions by combining significant SNPs within 10 Mb, and suggest these 

197 regions are likely to encompass the causative gene. 

198 As expected, the power to detect one or both causative loci depended on the location 

199 (autosomal vs. X-linked), dominance, and frequency of both 'causative' alleles (Figure 4, Figure 

200 4-figure supplement 4). For example, the mean percentage of simulations for which both loci 

201 were detected (SNP <10 Mb significant by permutation-based threshold) was six times higher 

202 (14.4%) for the X chromosome x autosomal dominant architecture compared to the autosomal- 

203 recessive x autosomal-recessive architecture (2.6%). The relationship between power and the 

204 proportion of affected individuals in the mapping population was complex. Interestingly, power 

p. 9 
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205 was high for some simulations with very few affected individuals. In these cases, the few 

206 individuals carrying the lower frequency sterility allele by chance also carried the sterility allele 

207 from the second locus, thus the average effect was not diminished by individuals carrying one 

208 but not both interacting sterility alleles. 

209 It is important to note that our empirical results suggest that the two-locus models used to 

210 simulate phenotypes are overly simplified. We predict that involvement of a sterility locus in 

211 multiple incompatibilities would reduce the influence of allele/genotype frequencies of any 

212 single partner locus on power. 

213 To estimate the false discovery rate from simulations, we classified significant SNPs not 

214 located on the same chromosome as either causative SNP as false positives. Choosing an 

215 appropriate distance threshold for false vs. true positives on the same chromosome was not 

216 obvious given the distribution of distances to causal SNPs (Figure 4-figure supplement 3). We 

217 classified significant SNPs <50 Mb from causative SNPs as true positives and excluded SNPs 

218 >50 Mb when calculating FDR. Using permutation-based significance thresholds, the median 

219 false positive rate was 0.014 (calculated for simulations with >10 SNPs within 50 Mb of either 

220 causative locus). These results suggest significant SNPs from the GWAS identified using this 

221 more stringent threshold are likely to be true positives. By contrast, the median false positive rate 

222 was 0.280 using the FDR<0.1 threshold, indicating this threshold is more permissive than 

223 predicted. Thus, there is a substantial chance that SNP associations with relative testis weight 

224 and expression PC 1 identified using this threshold are spurious and evidence is weak for GWAS 

225 regions comprising one SNP significantly associated with a single phenotype. 
226 

227 
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228 DISCUSSION 

229 Genetic mapping of testis weight and testis gene expression in hybrid zone mice implicated 

230 multiple autosomal and X-linked loci, and a complex set of interactions between loci. These 

23 1 results provide insight into the genetic architecture of a reproductive barrier between two 

232 incipient species in nature. 

233 Association mapping in natural hybrid populations. The potential to leverage 

234 recombination events from generations of intercrossing in hybrid zones to achieve high- 

235 resolution genetic mapping of quantitative traits has been recognized for decades (reviewed in 

236 Rieseberg and Buerkle 2002)1 . Until recently, collection of dense genotype datasets and large 

237 sample sizes has not been feasible in natural populations due to logistics and costs. This study 

238 demonstrates that loci and genetic interactions contributing to reproductive barrier traits can be 

239 identified in a GWAS with a modest sample size. Sample sizes approximating those used for 

240 human GWAS are not necessary if the prevalence and genetic architecture of the trait of interest 

241 are favorable. In general, epistasis makes genetic mapping more difficult. However, for hybrid 

242 defects, dependence of the phenotype on epistasis conversely may facilitate mapping. Despite 

243 substantial deleterious effects in hybrids, incompatibility alleles are not subject to negative 

244 selection within species and may be at high frequency or fixed within species. Hence, the 

245 prevalence of affected individuals in a hybrid zone for epistatic traits may be much higher than 

246 for deleterious traits in pure populations {e.g. disease in humans). 

247 Combining mapping of multiple sterility-related phenotypes substantially improved 

248 power to identify sterility loci. A few loci were identified for individual phenotypes using 

249 stringent significance thresholds. However, most loci identified using permissive thresholds 

250 showed significant associations with more than one phenotype. Spurious associations are 

p. 11 
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251 unlikely to be shared across phenotypes, thus evidence from multiple phenotypes provided 

252 confidence for contributions of 9 genomic regions to testis weight (on the X and 5 autosomes) 

253 and 50 genomic regions to expression PCI (on the X and 18 autosomes). 

254 The high resolution of mapping in the hybrid zone provides an advantage over laboratory 

255 crosses. For example, significant regions identified here (median = 2.1 Mb, regions with defined 

256 intervals) are much smaller than sterility QTL identified in F 2 s (White et al. 201 1). Many GWAS 

257 regions contain few enough genes that it will be possible to individually evaluate the potential 

258 role of each in future studies to identify causative genes. For example, 8/12 testis- weight regions 

259 and 28/50 expression PCI regions contain 10 or fewer protein-coding genes. We identified 

260 candidate genes with known roles in reproduction in several GWAS regions (Tables 1-2). 

261 However, for the majority of regions (8/12 relative testis weight, 33/50 expression PCI), there 

262 are no overlapping/nearby candidate genes. It is unlikely that these regions would be prioritized 

263 if contained in large QTL intervals. High resolution mapping is possible using mapping 

264 resources such as the collaborative cross (Aylor et al. 201 1) and heterogeneous stocks (Svenson 

265 et al. 2012), but these populations represent a small proportion of genetic diversity in house mice 

266 (Yang et al. 201 1) and hybrid incompatibility alleles may have been lost during strain 

267 production. 

268 Polymorphism of hybrid male sterility loci. Comparisons of different Fi crosses 

269 between strains of domesticus and musculus have shown that hybrid sterility phenotypes and loci 

270 depend on the geographic origins of parental strains (Britton-Davidian et al. 2005; Good et al. 

271 2008a), suggesting that most hybrid sterility alleles are segregating as polymorphisms within 

272 subspecies. Several of the loci identified in this study of hybrid zone mice are novel, providing 

273 additional evidence that sterility alleles are polymorphic within subspecies. However, a majority 
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274 of loci we identified in natural hybrids are concordant with previously identified sterility QTL 

275 (Tables 1-2, Figure 2). This similarity suggests there are common genetic factors underlying 

276 hybrid sterility in house mice, although there was no statistical support that genome-wide 

277 patterns of overlap with previous studies for testis weight or expression PC 1 were non-random 

278 (P > 0.05, 10,000 permutations). 

279 Prdm9, discovered by mapping Fi hybrid sterility, is the only characterized hybrid 

280 sterility gene in mice (Mihola et al. 2009). None of the GWAS regions identified here overlap 

281 Prdm9 (chromosome 17, 15.7 Mb). However, one expression PCI region (PC42) is ~4 Mb 

282 proximal to Prdm9. Reductions in PCI are observed in individuals that are heterozygous or 

283 homozygous for the domesticus allele at PC42. This pattern is partially consistent with sterility 

284 caused by Prdm9, which occurs in heterozygous individuals carrying sterile alleles from 

285 domesticus (Dzur-Gejdosova et al. 2012; Flachs et al. 2012). We did not find evidence for 

286 significant associations between SNPs near Prdm9 and testis weight; the nearest GWAS region 

287 (RTW09) is -41 Mb distal and low testis-weight is associated with the musculus allele. 

288 There is concordance between some of the genetic interactions between loci identified 

289 here and interactions identified by mapping sterility phenotypes and testis expression traits in an 

290 F 2 cross between musculus and domesticus (White et al. 201 1; Turner et al. 2014) (Figure 2 - 

291 figure supplement 1). Precise overlap between some GWAS regions and interaction regions from 

292 F 2 s identifies strong candidates for future studies to identify the causative mechanisms and genes 

293 underlying sterility loci. For example, an interaction between chromosome 12 and the central X 

294 chromosome (RTW1 1, PC49) identified for testis weight and expression PCI overlaps an 

295 interaction affecting testis expression in F 2 hybrids (Turner et al. 2014). The 4.3 Mb interval of 

296 overlap among chromosome 12 loci (RTW07, PC29, 32.38 - 41.43 Mb F 2 s) encompasses 12 

p. 13 
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297 protein-coding genes, including a gene with a knockout model showing low testis weight and 

298 sperm count (Arl4a) (Schurmann et al. 2002), and two genes with roles in regulating gene 

299 expression (Meox2, Etvl). 

300 We compared the positions of GWAS regions to 182 regions (163 autosomal, 19 X- 

301 linked) with evidence for epistasis based on a genome -wide analysis of genomic clines in a 

302 transect across the house mouse hybrid zone in Bavaria (Janousek et al. 2012), the same region 

303 where the progenitors of the mapping population were collected. Five testis-weight regions and 

304 18 expression-PC 1 regions overlap candidate regions from the hybrid zone genomic clines 

305 analysis (Tables 1-2), however the patterns of overlap were not statistically significant {P > 0.05, 

306 10,000 permutations). Future introgression analyses using high-density markers within and 

307 around GWAS regions may be useful in identifying causative genes and estimating the 

308 contributions of sterility alleles to reduced gene flow. 

309 Role of the X chromosome. Three GWAS regions associated with testis weight and five 

310 expression PCI regions are located on the X chromosome. The X-chromosomal regions surpass 

311 the stringent permutation-based significance threshold, and thus have strong statistical support. 

3 12 These results are consistent with evidence for an important role for the X in hybrid sterility from 

313 laboratory crosses between subspecies strains geographically diverse in origin (Guenet et al. 

314 1990; Elliott et al. 2001; Oka et al. 2004; Storchova et al. 2004; Oka et al. 2007; Good et al. 

315 2008a,b; Mihola et al. 2009; White et al. 2012) and evidence for greatly reduced gene flow of X- 

316 linked loci across the European hybrid zone (Tucker et al. 1992; Payseur et al. 2004; Macholan 

317 et al. 2007; Teeter et al. 2008; 2010). A disproportionately large contribution of the X 

318 chromosome is a common feature of reproductive isolation in many taxa, the so-called "large X 

319 effect" (Coyne and Orr 1989). 

p. 14 
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320 The musculus derived X chromosome has been implicated repeatedly in genetic studies 

321 of sterility in Fi and F 2 hybrids (reviewed in Good et al. 2008a; White et al. 201 1). By contrast, 

322 domesticus alleles were associated with the sterile pattern for most loci we identified on the X in 

323 hybrid zone mice (Tables 1-2). A testis expression-QTL mapping study performed in F 2 s showed 

324 that domesticus ancestry in the central/distal region of the X was associated with a sterile 

325 expression pattern (Turner et al. 2014). Differences between studies might reflect geographic 

326 variation in sterility alleles but identification of domesticus-sterilc X alleles only in generations 

327 beyond the Fi suggests interactions with recessive autosomal partner loci are essential. The 

328 importance of recessive sterility alleles was demonstrated by the discovery of multiple novel 

329 recessive loci in an F 2 mapping study (White et al. 201 1). Fi hybrids are essentially absent in 

330 nature (Teeter et al. 2008; Turner et al. 2012), because the hybrid zone is >30 km wide (Boursot 

331 et al. 1993) thus pure subspecies individuals rarely encounter each other. Consequently, 

332 recessive autosomal loci acting in later generations are essential in maintaining reproductive 

333 isolation at present and likely contributed to its establishment. 

334 Genetic architecture of hybrid sterility. Despite a growing list of sterility loci and 

335 genes identified in a variety of animal and plant taxa, there are few cases of Dobzhansky-Muller 

336 incompatibilities for which all partner loci are known (Phadnis 201 1). Hence, there remain many 

337 unanswered questions about the genetic architecture of hybrid defects. For example, how many 

338 incompatibilities contribute to reproductive barriers in the early stages of speciation? How many 

339 partner loci are involved in incompatibilities? Are these patterns consistent among taxa? 

340 The interactions contributing to sterility phenotypes we mapped in hybrid zone mice 

341 reveal several general features of the genetic architecture of hybrid sterility. Most sterility loci 

342 interact with more than one partner locus. This pattern is consistent with evidence from studies 

p. 15 
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343 mapping sterility in Fi musculus-domesticus hybrids (Dzur-Gejdosova et al. 2012) and mapping 

344 interactions affecting testis gene expression in F 2 hybrids (Turner et al. 2014). We did not have 

345 sufficient power to map interactions requiring three or more sterility alleles, but interactions 

346 between alleles from the same subspecies imply their existence. Loci causing male sterility in 

347 Drosophia pseudoobscura Bogota-USA hybrids also have multiple interaction partners; seven 

348 loci of varying effect size interact to cause sterility (Phadnis 201 1). These results suggest 

349 biological pathways/networks are often affected by multiple Dobzhansky-Muller interactions; a 

350 single pairwise interaction between incompatible alleles disrupts pathway function enough to 

351 cause a hybrid defect phenotype but when more incompatible alleles are present, the effects of 

352 multiple pairwise interactions are synergistic. Variation in the effect sizes of sterility loci might 

353 then reflect variation in the number of networks in which the gene is involved, and the 

354 connectedness/centrality of the gene within those networks. 

355 Characteristics of the incompatibility network are important for generating accurate 

356 models of the evolution of reproductive isolation. A "snowball effect" - faster-than-linear 

357 accumulation of incompatibilities caused by epistasis - is predicted on the basis of the 

358 Dobzhansky-Muller model (Orr 1995; Orr and Turelli 2001). Patterns of accumulation of hybrid 

359 incompatibilities in Drosophila and Solarium provide empirical support for the snowball 

360 hypothesis (Moyle and Nakazato 2010; Matute et al. 2010). Because most GWAS regions have 

361 many interaction partners, our results are not consistent with the assumption of the snowball 

362 model that incompatibilities are independent, suggesting network models of incompatibilities 

363 (Johnson and Porter 2000; Porter and Johnson 2002; Johnson and Porter 2007; Palmer and 

364 Feldman 2009) may be more accurate for understanding the evolution of reproductive barriers in 

365 house mice. 
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366 Involvement of hybrid sterility loci in interactions with multiple partner loci also has 

367 important implications for understanding the maintenance of the hybrid zone. Because 

368 deleterious effects of a sterility allele are not dependent on a single partner allele, the marginal 

369 effect of each locus and thus visibility to selection is less sensitive to the allele frequencies at any 

370 single partner locus in the population. 

371 Identifying and functionally characterizing incompatibility genes is an important goal in 

372 understanding speciation, but is unrealistic in most non-model organisms. By contrast, mapping 

373 reproductive isolation traits in natural populations to identify the number and location of loci and 

374 interactions is feasible. General features of the genetic architecture of hybrid sterility - the 

375 number of incompatibilities and number and effect size of interacting loci - are arguably more 

376 likely to be shared among organisms than specific hybrid sterility genes. Comparison of these 

377 features among taxa may reveal commonalities of the speciation process. 
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378 MATERIALS AND METHODS 

379 Mapping population. The mapping population includes first-generation lab-bred male 

380 offspring of mice captured in the hybrid zone (Bavaria) in 2008 (Turner et al. 2012) (Figure 1 - 

381 figure supplement 1). We included 185 mice generated from 63 mating pairs involving 37 

382 unrelated females and 35 unrelated males. Many dams and sires were used in multiple mating 

383 pairs, thus our mapping population includes full siblings, half siblings and unrelated individuals. 

384 Most mating pairs (53 pairs, 149 offspring) were set up with parents originating from the same or 

385 nearby trapping locations. Eleven pairs (36 offspring) include dams and sires originating from 

386 more distant trapping locations; phenotypes of these offspring were not reported in (Turner et al. 

387 2012). 

388 Phenotyping. Males were housed individually after weaning (28 days) to prevent effects 

389 of dominance interactions on fertility. We measured combined testis weight and body weight 

390 immediately after mice were sacrificed at 9-12 weeks. We calculated relative testis weight (testis 

391 weight/body weight) to account for a significant association between testis weight and body 

392 weight (Pearson's correlation = 0.29, P = 4.9xl0" 5 ). 

393 We classify individuals with relative testis weight below the range observed in pure 

394 subspecies as showing evidence for sterility (Turner et al. 2012). To confirm that this is an 

395 appropriate threshold for inferring hybrid defects, we compared this value to relative testis 

396 weights reported previously for offspring from intraspecific and interspecific crosses (Good et al. 

397 2008a). The pure subspecies minimum we observed is substantially lower (>2 standard 

398 deviations) than means for males from intraspecific crosses (converted from single relative testis 

399 weight: musculus ?WK x musculus CZECH - mean = 10.2, standard deviation = 1 .2; domesticus LEWES 

400 x domesticus™ 3 - mean = 1 1.0, standard deviation = 1.0) and comparable to (within 1 standard 
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401 deviation) values observed in Fi hybrids from 4/7 interspecific crosses that showed significant 

402 reductions (mean plus one standard deviation 4.6 - 9.2mg/g). 

403 Testis gene expression. We measured gene expression in testes of 179 out of the 185 

404 males from the mapping population. Freshly dissected testes were stored in RNAlater (Qiagen) at 

405 4° overnight, then transferred to -20° until processed. We extracted RNA from 15-20 mg whole 

406 testis using Qiagen RNeasy kits, and a Qiagen Tissue Lyser for the homogenization step. We 

407 verified quality of RNA samples (RIN > 8) using RNA 6000 Nano kits (Agilent) on a 2100 

408 Bioanalyzer (Agilent). 

409 We used Whole Mouse Genome Microarrays (Agilent) to measure genome-wide 

410 expression. This array contains 43,379 probes surveying 22,210 transcripts from 21,326 genes. 

411 We labeled, amplified, and hybridized samples to arrays using single-color Quick- Amp Labeling 

412 Kits (Agilent), according to manufacturer protocols. We verified the yield (>2 |ig) and specific 

413 activity (>9.0 pmol Cy3/|lg cRNA) of labeling reactions using a NanoDrop ND-1000 UV-VIS 

414 Spectrophotometer (NanoDrop, Wilimington, DE, USA). We scanned arrays using a High 

415 Resolution Microarray Scanner (Agilent) and processed raw images using Feature Extraction 

416 Software (Agilent). Quality control procedures for arrays included visual inspection of raw 

417 images and the distribution of non-uniformity outliers to identify large spatial artifacts (e.g. 

418 caused by buffer leakage or dust particles) and quality control metrics from Feature Extraction 

419 protocol GEl_QCMT_Dec08. 

420 We mapped the 41,174 non-control probe sequences from the Whole Mouse Genome 

421 Microarray to the mouse reference genome (NCBI37, downloaded March 201 1) using BLAT 

422 ((Kent 2002); minScore = 55, default settings for all other options). Probes with multiple perfect 

423 matches, more than nine imperfect matches, matches to non-coding/intergenic regions only, or 
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424 matches to more than one gene were excluded. A total of 36,323 probes (covering 19,742 Entrez 

425 Genes) were retained. 

426 We preformed preprocessing of microarray data using the R package Agi4x44PreProcess 

427 (Lopez-Romero 2009). We used the background signal computed in Feature Extraction, which 

428 incorporates a local background measurement and a spatial de-trending surface value. We used 

429 the "half setting in Agi4x44PreProcess, which sets intensities below 0.5 to 0.5 following 

430 background subtraction, and adds an offset value of 50. Flags from Feature Extraction were used 

43 1 to filter probes during preprocessing (wellaboveBG=TRUE, isfound=TRUE, 

432 wellaboveNEG=TRUE). We retained probes with signal above background for at least 10% of 

433 samples. We used quantile normalization to normalize signal between arrays. Expression data 

434 were deposited in Gene Expression Omnibus as project GSE61417. 

435 To identify major axes of variation in testis expression, we performed a principal 

436 components analysis using prcomp in R (R Development Core Team 2010) with scaled 

437 variables. 

438 Genotyping. We extracted DNA from liver, spleen, or ear samples using salt extraction 

439 or DNeasy kits (Qiagen, Hilden, Germany). Males from the mapping population were geno typed 

440 using Mouse Diversity Genotyping Arrays (Affymetrix, Santa Clara, CA) by Atlas Biolabs 

441 (Berlin, Germany). 

442 We called genotypes at 584,729 SNPs using apt-probeset-genotype (Affymetrix) and 

443 standard settings. We used the MouseDivGeno algorithm to identify variable intensity 

444 oligonucleotides (VINOs) (Yang et al. 201 1); 53,148 VINOs were removed from the dataset. In 

445 addition, we removed 18,120 SNPs with heterozygosity > 0.9 in any population because these 

446 SNPs likely represent additional VINOs. We performed additional filtering steps on SNPs 
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447 included in the dataset used for mapping. We only included SNPs with a minor allele frequency 

448 >5% in the mapping population. SNPs without a genome position or with missing data for >15% 

449 of the individuals in the mapping population or pure subspecies reference panel were removed. 

450 We pruned the dataset based on linkage disequilibrium (LD) to reduce the number of tests 

451 performed. LD pruning was performed in PLINK (Purcell et al. 2007; Purcell n.d.) using a 

452 sliding window approach (30 SNPs window size, 5 SNPs step size) and a VIF threshold of 1 x 

453 10" 6 (VIF = 1/(1 -R 2 ) where R 2 is the multiple correlation coefficient for a SNP regressed on all 

454 other SNPs simultaneously). This procedure essentially removed SNPs in perfect LD. These 

455 filtering steps yielded 156,204 SNPs. 

456 Ancestry inference. To identify ancestry-informative SNPs, we compared genotypes 

457 from 21 pure M. m. domesticus individuals (11 from Massif Central, France and 10 from 

458 Cologne/Bonn, Germany) and 22 M. m. musculus individuals (11 from Namest nad Oslavou, 

459 Czech Republic and 1 1 from Almaty, Kazakhstan) (Staubach et al. 2012). 

460 We used Structure (Pritchard et al. 2000; Falush et al. 2003) to graphically represent the 

461 genetic composition of our mapping population (Figure 1 -figure supplement 1). We included one 

462 diagnostic SNP per 20 cM, 3-5 markers/chromosome totaling 60 SNPs genome wide. We used 

463 the 'admix' model in Structure and assumed two ancestral populations. 

464 Association mapping. To identify genomic regions significantly associated with relative 

465 testis weight and testis gene expression, we used a mixed model approach to test for single SNP 

466 associations. Admixture mapping - often applied in studies using samples with genetic ancestry 

467 from two distinct populations - was not appropriate for this study because it was not possible to 

468 account for relatedness among individuals in the mapping population (Buerkle and Lexer 2008) 

469 (Winkler et al. 2010). 
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470 We performed association mapping using GEMMA (Zhou and Stephens 2012), which 

471 fits a univariate mixed model, incorporating an n x n relatedness (identity-by-state) matrix as a 

472 random effect to correct for genetic structure in the mapping population. We estimated 

473 relatedness among the individuals in the mapping population in GEMMA using all markers and 

474 the -gk 1 option, which generates a centered relatedness matrix. For each single SNP association 

475 test we recorded the Wald test P value. Phenotypes tested include relative testis weight (testis 

476 weight/body weight, RTW), testis expression principal component 1 (PCI, 14.6% variance, 

477 associated with fertility, Figure 1 - figure supplement 2), and normal quantile ranks of gene 

478 expression values for individual transcripts. Neither RTW nor expression PCI were significantly 

479 correlated with age at phenotyping (RTW - cor= -0.02, P=0.72; PCI - cor= 0.01, P=0.90), thus 

480 we did not include age in the model. SNP data, phenotypic data and kinship matrix to run 

481 GEMMA area available through Dryad at: doi:l 0.506 l/dryad.2br40. 

482 To account for multiple testing, we used a significance threshold based on the 10% false 

483 discovery rate (Benjamini and Hochberg 1995), equivalent to P = 3.49 x 10" 5 for RTW and P = 

484 2.86 x 1 0" 4 for expression PC 1 . We determined more stringent significance thresholds by 

485 permutation. We randomized phenotypes among individuals 10,000 times, recording the lowest 

486 P value on the X and autosomes for each permutation. Thresholds set to the 5th percentile across 

487 permutations for RTW were 5.73 x 10" 7 (autosomes) and 5.83 x 10" 5 (X chromosome); thresholds 

488 for expression PCI were 1.66 x 10" 8 (autosomes) and 1.01 x 10" 5 (Xchromosome). 

489 To estimate the genomic interval represented by each significant LD-filtered SNP, we 

490 report significant regions defined by the most distant flanking SNPs in the full dataset showing r 

491 > 0.9 (genotypic LD, measured in PLINK) with each significant SNP. We combined significant 

492 regions < 10 Mb apart into a single region. 

p. 22 



Downloaded from http://biorxiv.org/on September 18, 2014 

Turner & Harr 

493 Testing for genetic interactions. Identifying genetic interactions using GWAS is 

494 computationally and statistically challenging. To improve power, we reduced the number of tests 

495 performed by testing for interactions only among significant SNPs (FDR < 0. 1) identified using 

496 GEMMA. We tested all pairs of significant SNPs located on different chromosomes for each 

497 phenotype (692 pairs RTW, 82,428 pairs expression PCI). To account for relatedness among 

498 individuals we used a mixed model approach, similar to the model implemented in GEMMA. 

499 We used the Imekin function from the coxme R package (Therneau 2012) to fit linear mixed 

500 models including the identity-by-state kinship matrix as a random covariate. We report 

501 interactions as significant for SNP pairs with PO.05 and FDR<0. 1 for interaction terms (RTW: 

502 FDRO.l ~ PO.02; expression PCI: FDRO.09 ~ PO.05). 

503 Mapping Simulations. We performed simulations to evaluate the performance of our 

504 mapping approach under varying genetic architectures and allele frequencies. We simulated 

505 phenotypes using several genetic models of two-locus epistasis and parameters based on the 

506 empirical distribution of relative testis weight. The simulation procedure is illustrated in Figure 4 

507 - figure supplement 1 . To preserve genetic structure, we simulated phenotypes using two-locus 

508 genotypes from the SNP dataset. 

509 We tested 100 autosomal-autosomal SNP pairs (SNPs on different chromosomes) and 

510 100 X- Autosomal pairs (50 with domesticus X-linked sterile alleles and 50 with musculus X- 

511 linked sterile alleles). The criteria used for choosing 'causative' SNPs were a minor allele 

512 frequency > 0.05 in the mapping population and fixed in at least one pure subspecies. The 

513 'sterile' allele could be polymorphic or fixed within subspecies but the alternate 'non-sterile' 

514 allele had to be fixed within the other subspecies -e.g. domesticus sterile alleles have frequencies 

515 0.05-1 .0 in the domesticus reference populations from France and Germany and the alternate 
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516 allele at those SNPs are fixed in mus cuius samples from the Czech Republic and Kazakhstan. For 

517 each pair, the 'causative' SNPs were randomly selected from all SNPs meeting those criteria 

518 (144,506 possible domesticus sterile, 124,390 possible musculus sterile). 

519 For each SNP pair, we modeled all possible combinations of recessive, additive, and 

520 dominant autosomal sterile alleles. For each model type, we assigned mean Z scores for each 

521 possible two-locus genotype (Figure 4 - figure supplement 2). The magnitude of the most severe 

522 phenotype (-2.3 standard deviations) is based on observed relative testis weights in the most 

523 severely affected males. The mean Z score for heterozygotes in additive models was -1.15. Mean 

524 Z scores for non-sterile genotypes in the models were randomly drawn from a uniform 

525 distribution between -0.5 and 0.5. 

526 For each SNP pair/architecture, 100 data sets were generated by drawing phenotypes (Z 

527 scores) for each individual from a normal distribution with the appropriate two-locus mean and 

528 standard deviation = 0.75. The standard deviation value, equivalent to 2.98 mg/g, was chosen on 

529 the basis of standard deviations in pure subspecies samples from the mapping population 

530 (domesticus = 2.13, musculus = 3.65; (Turner et al. 2012)). This value is higher than standard 

531 deviations in intraspecific Fl males (domesticus LEWES x domesticus WSB = 1.2, musculus™*^ x 

532 musculus CZECR = 1 .0; (Good et al. 2008a)), suggesting estimates of mapping power may be 

533 conservative. 

534 In total, 90,000 simulations were performed, (9 architectures x 100 SNP pairs x 100 data 

535 sets). We identified significant SNPs for each data set using GEMMA, as described above for the 

536 empirical data. Testing all pairwise interactions between significant SNPs for each simulated 

537 dataset was not feasible computationally. For each SNP pair and architecture, we tested all 

538 pairwise interactions for one randomly chosen replicate with 50 - 500 significant SNPs, a range 
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539 encompassing the observed number of significant associations for relative testis weight and 

540 expression PCI. Of the 900 SNP pair/architecture combinations, 889 had at least one replicate 

541 with results in this range. 

542 Significance of overlap between candidate sterility loci. We used permutations to test 

543 for non-random co-localization of candidate sterility loci from this study and previous QTL and 

544 hybrid zone studies. The locations of significant GWAS regions for relative testis weight and 

545 expression PCI were randomized 10,000 times using BEDTools (Quinlan and Hall 2010). To 

546 assess overlap between significant regions for the two phenotypes, we counted the number of 

547 RTW regions overlapping PCI regions (and vice versa) for each permutation. To test for overlap 

548 between GWAS identified regions and previously reported candidate regions for related 

549 phenotypes, we counted the number of permuted regions overlapping the positions of the 

550 published regions (fixed) for each replicate. GWAS regions for both phenotypes were compared 

551 to genomic regions with evidence for epistasis and reduced introgression in the Bavarian transect 

552 of the hybrid zone (Janousek et al. 2012). In addition, RTW regions were compared to testis 

553 weight QTL from mapping studies in F 2 and backcross hybrids from crosses between subspecies 

554 (Storchova et al. 2004; Good et al. 2008b; White et al. 201 1; Dzur-Gejdosova et al. 2012) and 

555 expression PCI regions were compared to trans eQTL hotspots identified in F 2 hybrids (Turner 

556 etal. 2014). 

557 Gene annotation. We used ENSEMBL (version 66, February 2012) Biomart to 

558 download gene annotations for genomic regions significantly associated with relative testis 

559 weight. We identified candidate genes in significant regions with roles in male reproduction 

560 using reviews of male fertility (Matzuk and Lamb 2008), manual searches, MouseMine searches 

561 for terms related to male fertility ( http : / / www . mousemine .org/ ) and gene ontology (GO) terms 
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562 related to male reproduction or gene regulation (plus children): meiosis GO:0007126; DNA 

563 methylation GO:0006306; regulation of gene expression GO:0010468; transcription 

564 GO:0006351; spermatogenesis GO:0007283; male gamete generation GO:0048232; gamete 

565 generation GO:0007276; meiotic cell cycle GO:0051321. Many genes with roles in reproduction 

566 reported in publications were not annotated with related GO terms, highlighting the limitations of 

567 gene ontology. Moreover, genes causing sterility might not have functions obviously related to 

568 reproduction. 
569 
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809 Figure legends: 

810 Figure 1. Single SNPs associated with (A) relative testis weight, (B) testis expression principal 

811 component 1, and (C) expression of transcripts located on other chromosomes (trans). Dashed 

812 lines indicate significance thresholds based on: permutations for autosomes (labeled 5% perm 

813 A), permutations for X chromosome (labeled 5% perm X), false discovery rate < 0.1 (labeled 

814 10% FDR), and 95 th percentile of significant transcript association counts across SNPs (labeled 

815 95%). 
816 

817 Figure 1 - Figure supplement 1. Mapping population. (A) Location of sampling area (black 

818 box) in European house mouse hybrid zone. (B) Sampling locations for parents of mice in the 

819 mapping population. (C) Structure analysis of mapping population. Individuals (vertical bands) 

820 are arranged by geographic origin and average percentage alleles from Mus musculus musculus. 
821 

822 Figure 1 - Figure supplement 2. Principal components analysis of genome-wide gene 

823 expression in testis. (A) Plot of principal component 1 (PCI) vs. PC2 scores. Individuals with 

824 relative testis weight and/or sperm count below the pure subspecies range are indicated in blue 

825 ("low fertility"). Individuals with relative testis weight and sperm count within one standard 

826 deviation of the mean in pure subspecies individuals are indicated in red ("fertile range"). (B) 

827 Plot of relative testis weight vs. PCI score. Correlation coefficient (Pearson's) and P value are 

828 indicated. "). (C) Plot of hybrid index (% musculus alleles on autosomes) vs. PC2 score. 

829 Correlation coefficient (Pearson's) and P value are indicated. 
830 
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831 Figure 1 - Source data 1. SNPs significantly associated with relative testis weight and/or testis 

832 expression PCI (excel file). 

833 

834 Figure 2. Significant GWAS regions and interactions associated with (A) relative testis weight 

835 and (B) testis expression principal component 1 in hybrid zone mice. In (A), orange and yellow 

836 boxes in outer rings (outside grey line) indicate quantitative trait loci (QTL) identified for testis 

837 weight and other sterility phenotypes in previous studies (see Table 1 for details). Green boxes 

838 indicate significant GWAS regions for relative testis weight. Green lines represent significant 

839 genetic interactions between regions; shade and line weight indicate the number of significant 

840 pairwise interactions between SNPs for each region pair. In (B), orange boxes in outer rings 

841 indicate QTL for testis-related phenotypes (testis weight and seminiferous tubule area) identified 

842 in previous studies, yellow boxes indicate QTL for other sterility phenotypes and red boxes 

843 indicate trans eQTL hotspots (see Table 2 for details). Green boxes indicate significant GWAS 

844 regions for relative testis weight. Purple boxes indicate significant GWAS regions for testis 

845 expression PCI. Lines represent significant genetic interactions between regions; color and line 

846 weight - as specified in legend - indicate the number of significant pairwise interactions between 

847 SNPs for each region pair. Plot generated using circos (Krzywinski et al. 2009). 
848 

849 Figure 2 - Figure supplement 1. Genetic interactions associated with hybrid sterility hybrid 

850 zone mice and in F2 hybrids. Orange boxes in outer rings indicate QTL for testis-related 

851 phenotypes (testis weight and seminiferous tubule area) identified in previous studies, yellow 

852 boxes indicate QTL for other sterility phenotypes and red boxes indicate trans eQTL hotspots 

853 (see Table 2 for details). Green boxes indicate significant GWAS regions for relative testis 
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854 weight. Purple boxes indicate significant GWAS regions for testis expression PCI. Lines 

855 represent significant genetic interactions identified in hybrid zone mice for relative testis weight 

856 (in green) and expression PCI (in purple) which are concordant with genetic interactions 

857 identified by mapping expression traits in F2 hybrids (Turner et al. 2014). Plot generated using 

858 circos (Krzywinski et al. 2009). 
859 

860 Figure 2 - Source data 1. Significant genetic interactions (SNP pairs) for relative testis weight 

861 (excel file). 

862 

863 Figure 2 - Source data 2. Significant genetic interactions (SNP pairs) for testis expression PCI 

864 (excel file). 

865 

866 Figure 3. Phenotypic effects of testis-weight loci and interactions. Histograms showing 

867 maximum deviations from the population mean for (A) single SNPs and (B) two-locus 

868 interactions. Dashed vertical lines indicate minimum values observed in pure subspecies males. 

869 (C) Examples of phenotypic means by two-locus genotype for autosomal-autosomal and X- 

870 autosomal interactions. Genotypes are indicated by one letter for each locus: D - homozygous 

871 for the domesticus allele, H - heterozygous, M - homozygous musculus. 
872 

873 Figure 4. Mapping power in simulations. Each panel results from a single genetic architecture 

874 model for (A) 100 autosomal-autosomal SNP pairs and (B) 100 X-autosomal SNP pairs. Each 

875 point represents the percentage of data sets generated from a single SNP pair in which locus 1 

876 (domesticus sterile allele; green), locus 2 (musculus sterile allele; purple), or both loci (orange) 
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877 were identified by association mapping (>1 SNP significant by permutation based threshold 

878 within 10 Mb of 'causal' SNP). The x axis indicates the percentage of individuals with partial or 

879 full sterility phenotypes. Curves were fit using 2 nd order polynomials. In (A), locus 1 indicates 

880 the SNPs with musculus alleles sterile and locus 2 indicates the SNPs with domes ticus alleles 

881 sterile. In (B), locus 1 is the X-linked SNP and locus 2 is the autosomal SNP. 
882 

883 Figure 4 - Figure Supplement 1. Mapping simulation methods. Schematics of (A) choice of 

884 'causal' SNP pairs from the genotype data, (B) phenotype distributions for simulations, (C) 

885 generation of simulated phenotype data sets, (D) association mapping. In (B), histogram shows 

886 the empirical distribution of relative testis weight in the mapping population, in standard 

887 deviation units. 
888 

889 Figure 4 - Figure Supplement 2. Table: Z scores for simulation models. 

890 

891 Figure 4 - Figure Supplement 3. Distances of significant SNPs to causal SNP in simulations. 

892 Distributions are shown at two scales for autosomal and X-linked loci. 

893 

894 Figure 4 - Figure Supplement 4. Table: Results of mapping simulations. 
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895 Table 1. Genomic regions significantly associated with relative testis weight. 
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897 "Significant intervals were defined by positions of the most proximal and distal SNPs with LD > 0.9 to a significant SNP. 

898 3 The number of SNPs significant at FDR <0. 1 is reported; number of significant SNPs significant with <0.05 P value in permutations is in parentheses. 

899 4 Number of significant SNPs enriched for associations with transcripts expressed on another chromosome (P <0.05; FDR<0. 1; >30 transcripts). 

900 5 Number of regions with significant interactions. 

901 Overlapping regions significant for expression PCI (see Table 2). 

902 7 Sterility QTL overlapping or within 10 Mb from A (White et al. 201 1), B (Dzur-Gejdosova et al. 2012), c (Turner et al. 2014), D (Good et al. 2008b). Abbreviations 

903 for phenotypes: ASH: abnormal sperm head morphology, TW: testis weight, SC: sperm count, shPCl: sperm head shape PCI, eQTLHS: trans eQTL hotspot 
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904 tubule area, FERT: fertility, PBT: proximal bent sperm tail, HT: headless/tailless sperm, DBT: distal bent sperm tail, TAS: total abnormal sperm. BHZ: 

905 overlapping candidate regions with evidence from epistasis in the Bavarian hybrid zone transect (Janousek et al. 2012). 

906 8 Sterile allele inferred on the basis of frequency of a majority of significant SNPs in pure subspecies samples: D - domesticus; M - musculus; Lower-case 

907 indicates F ST < 0.7 between pure subspecies; * indicates overlapping PCI region is D sterile; U - nondiagnostic SNP and/or no majority allele. 

908 9 Number of genes (protein-coding) overlapping region. 

909 10 Genes with roles in male reproduction on the basis of E male reproduction gene ontology terms (see Methods) or phenotypes of knockout models reported in 

910 F (Matzuk and Lamb 2008) or G MGI database. 
911 

912 Table 1 - Source data 1. Protein-coding genes in significant relative testis weight regions. 



913 Table 2. 
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914 'Significant SNPs <10 Mb apart were combined into regions. 

915 Significant intervals were defined by positions of the most proximal and distal SNPs with LD > 0.9 to a significant SNP. 

916 3 The number of SNPs significant at FDR <0.1 is reported; number of significant SNPs significant with <0.05 P value in permutations is in parentheses. 

917 4 Number of significant SNPs enriched for associations with transcripts expressed on another chromosome (P <0.05; FDR<0. 1; >30 transcripts). 

918 5 Number of regions with significant interactions. 

919 Overlapping regions significant for expression PCI (see Table 2). 

920 7 Sterility QTL overlapping or within 10 Mb from A (White et al. 201 1), B (Dzur-Gejdosova et al. 2012), c (Turner et al. 2014), D (Good et al. 2008b), E (Storchova et 

921 al. 2004). Abbreviations for phenotypes: ASH: abnormal sperm head morphology, TW: testis weight, SC: sperm count, shPCl: sperm head shape PCI, eQTLHS: 

922 trans eQTL hotspot, STA: seminiferous tubule area, FERT: fertility, PBT: proximal bent sperm tail, HT: headless/tailless sperm, DBT: distal bent sperm tail, 

923 TAS: total abnormal sperm, OFF: number of offspring. BHZ: overlapping candidate regions with evidence from epistasis in the Bavarian hybrid zone transect 

924 (Janousek et al. 2012). 

925 8 Sterile allele inferred on the basis of frequency of a majority of significant SNPs in pure subspecies samples: D - domesticus; M - musculus; Lower-case 

926 indicates F ST < 0.7 between pure subspecies; * indicates overlapping PCI region is D sterile; U - nondiagnostic SNP and/or no majority allele; Dh - two SNPs 

927 with domesticus sterile alleles, one SNP heterozygous genotype shows sterile pattern; Md - majority musculus sterile alleles but some SNPs diagnostic 

928 domesticus sterile alleles. 

929 'Number of genes (protein-coding) overlapping region. 

930 10 Genes with roles in male reproduction on the basis of F male reproduction gene ontology terms (see Methods) or phenotypes of knockout models reported in 

93 1 G (Matzuk and Lamb 2008] or H MGI database. 

932 Table 2 - Source data 1. Protein-coding genes in significant testis expression PCI regions. 

933 
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